nanopyx.core.transform.sr_temporal_correlations

View Source

  1import numpy as np
  2
  3
  4def calculate_SRRF_temporal_correlations(im: np.array, order: int = 1, do_integrate_lag_times: bool = 0): #these are for SRRF
  5
  6    im = np.array(im, dtype='float32')
  7    assert im.ndim == 3 and order <= 4
  8    
  9    if order == 0: #TODO: check order number in NanoJ
 10        out_array = np.amax(im, axis=0)
 11    
 12    elif order == 1:
 13        out_array = np.mean(im, axis=0)
 14    
 15    elif order == -1:
 16        out_array = calculate_pairwise_product_sum(im)
 17
 18    else: # order = 2 or order = 3 or order = 4
 19        out_array = calculate_acrf_(im, order, do_integrate_lag_times)
 20    
 21    return out_array
 22
 23
 24def calculate_eSRRF_temporal_correlations(im: np.array, correlation: str):
 25
 26    im = np.array(im, dtype='float32')
 27
 28    if correlation == "AVG":
 29        out_array = np.mean(im, axis=0)
 30
 31    elif correlation == "VAR":
 32        out_array = np.var(im, axis=0)
 33
 34    elif correlation == "TAC2":
 35        out_array = calculate_tac2(im)
 36
 37    else:
 38        raise ValueError(f"Type of correlation must be AVG, VAR or TAC2")
 39    
 40    return out_array
 41
 42
 43    
 44def calculate_pairwise_product_sum(rad_array):
 45    n_time_points, height_m, width_m = rad_array.shape
 46
 47    out_array = np.zeros((height_m, width_m), dtype=np.float32)
 48    #max_array = np.max(rad_array, axis=0)
 49    counter = 0
 50    pps = 0
 51
 52    for t0 in range(n_time_points):
 53        r0 = np.maximum(rad_array[t0],0)
 54        if np.any(r0) > 0:
 55            for t1 in range(t0, n_time_points):
 56                r1 = np.maximum(rad_array[t1],0)
 57                pps += r0 * r1
 58                counter += 1
 59        else:
 60            counter += n_time_points - t0
 61    pps = pps/ max(counter,1)
 62    out_array = pps
 63
 64    return out_array
 65
 66
 67def calculate_acrf_(rad_array, order, do_integrate_lag_times):
 68    im = rad_array.copy()
 69    n_time_points, height_m, width_m = im.shape
 70    mean = np.mean(im, axis=0)
 71
 72    out_array = np.zeros((height_m, width_m), dtype=np.float32)
 73
 74    abcd = np.zeros((height_m, width_m), dtype=np.float32)
 75    abc = np.zeros((height_m, width_m), dtype=np.float32)
 76    ab = np.zeros((height_m, width_m), dtype=np.float32)
 77    cd = np.zeros((height_m, width_m), dtype=np.float32)
 78    ac = np.zeros((height_m, width_m), dtype=np.float32)
 79    bd = np.zeros((height_m, width_m), dtype=np.float32)
 80    ad = np.zeros((height_m, width_m), dtype=np.float32)
 81    bc = np.zeros((height_m, width_m), dtype=np.float32)
 82
 83    if do_integrate_lag_times != 1:
 84        t = 0
 85        while (t < n_time_points - order):
 86            ab = ab + (im[t] - mean) * (im[t+1] - mean)
 87            if order == 3:
 88                abc = abc + (im[t] - mean) * (im[t+1] - mean) + (im[t+2] - mean)
 89            if order == 4:
 90                a = im[t] - mean
 91                b = im[t+1] - mean
 92                c = im[t+2] - mean
 93                d = im[t+3] - mean
 94                abcd = abcd + np.multiply(np.multiply(np.multiply(a,b),c),d)
 95                cd = cd + np.multiply(c,d)
 96                ac = ac + np.multiply(a,c)
 97                bd = bd + np.multiply(b,d)
 98                ad = ad + np.multiply(a,d)
 99                bc = bc + np.multiply(b,c)
100            t = t + 1
101        if order == 3:
102            out_array = np.absolute(abc) / n_time_points
103        elif order == 4:
104            out_array = np.absolute(abcd - ab * cd - ac * bd - ad * bc) / n_time_points
105        else:
106            out_array = np.absolute(ab) / n_time_points
107    
108    else:
109        n_binned_time_points = n_time_points
110        tbin = 0
111        while n_binned_time_points > order:
112            t = 0
113            ab = np.zeros((height_m, width_m), dtype=np.float32)
114            while (t < n_binned_time_points - order):
115                tbin = t * order
116                ab = ab + (im[t] - mean) * (im[t+1] - mean)
117                if order == 3:
118                    abc = abc + (im[t] - mean) * (im[t+1] - mean) + (im[t+2] - mean)
119                if order == 4:
120                    a = im[t] - mean
121                    b = im[t+1] - mean
122                    c = im[t+2] - mean
123                    d = im[t+3] - mean
124                    abcd = abcd + np.multiply(np.multiply(np.multiply(a,b),c),d)
125                    cd = cd + np.multiply(c,d)
126                    ac = ac + np.multiply(a,c)
127                    bd = bd + np.multiply(b,d)
128                    ad = ad + np.multiply(a,d)
129                    bc = bc + np.multiply(b,c)
130                    
131                im[t] = np.zeros((height_m, width_m), dtype=np.float32)
132
133                if tbin < n_binned_time_points:
134                    for _t in range(order-1):
135                        im[t] = im[t] + np.divide(im[tbin + _t], order)
136                t = t + 1
137            if order == 3:
138                out_array = np.absolute(abc) / n_binned_time_points
139            elif order == 4:
140                out_array = np.absolute(abcd - ab * cd - ac * bd - ad * bc) / n_binned_time_points
141            else:
142                out_array = np.absolute(ab) / n_binned_time_points
143
144            n_binned_time_points = n_binned_time_points / order
145
146    return out_array
147
148
149def calculate_tac2(rad_array):
150    mean = np.mean(rad_array, axis=0)  
151    centered = rad_array - mean  # center data around the mean
152    nlag = 1  # number of lags to compute TAC2 for
153    out_array = np.mean(centered[:-nlag] * centered[nlag:], axis=0)
154
155    return out_array
156
157
158
159
160# class TemporalCorrelation:
161#     accepted_correlation_types = ('mean', 'var', 'tac2', 'mip', 'pair_wise_product_sum')
162
163#     def __init__(self, correlation_type: str):
164#         """
165#         Perform a temporal correlation analysis of an image (with shape time, height, width) #TODO: discuss what shape should be the input
166#         :param correlation_type: desired type of interpolation ("mean", "var" or "tac2")
167#         """
168#         self.correlation_type = correlation_type
169
170
171#     def calculate_tc(self, im: np.array):
172
173#         out_array = np.empty((im.shape[0]))
174
175#         # assert isinstance(im, np.ndarray)
176
177#         if im.dtype != np.float32:
178#             im = np.array(im,dtype='float32')
179
180#         assert im.ndim == 3 #TODO: discuss if we should consider (t,c,z,r,c) here
181
182#         if self.correlation_type == "mean":
183#             out_array = np.mean(im, axis=0)
184        
185#         elif self.correlation_type == "var":
186#             out_array = np.var(im, axis=0)
187        
188#         elif self.correlation_type == "tac2": # second order autocorrelation function
189#             mean = np.mean(im, axis=0)  
190#             centered = im - mean  # center data around the mean
191#             nlag = 1  # number of lags to compute TAC2 for
192#             out_array = np.mean(centered[:-nlag] * centered[nlag:], axis=0)
193
194#         elif self.correlation_type == "mip": # maximum intensity projection
195#              out_array = np.amax(im, axis=0) #calculate_mip(im, doIntensityWeighting)
196
197#         elif self.correlation_type == "pair_wise_product_sum":
198#             out_array = calculate_pairwise_product_sum_(im)
199        
200#         else:
201#              raise ValueError(f"Type of correlation must be one of {self.accepted_correlation_types}")
202
203#         return out_array
204    
205    
206    # def calculate_mip(im: np.array, doIntensityWeighting: bool = True):
207    #     """
208    #     Calculate Maximum Intensity Projection of an input image
209    #     """
210    #     #nb_time_points, height, width = im.shape
211    #     out = np.empty((im.shape[0]))
212
213    #     if doIntensityWeighting == True:
214    #         out = 
215    #     else:
216    #         out = np.amax(im, axis=0)
217    #     return out
218
219    # def calculate_pairwise_product_sum(im: np.array):
220    #     """
221    #     Calculate Pair-Wise Product Sum of a 3D input image
222    #     """
223
224
225
226    # def calculate_pairwise_product_sum(self, rad_array: np.array):
227    #     n_time_points, height_m, width_m = rad_array.shape
228
229    #     SRRF_array = np.zeros((height_m, width_m), dtype=np.float32)
230
231    #     r0 = np.maximum(rad_array[:, :, None, :], 0) 
232    #     r1 = np.maximum(rad_array[:, :, :, None], 0) 
233
234    #     pps = np.sum(r0 * r1, axis=(2, 3)) / np.maximum((n_time_points * (n_time_points + 1)) / 2, 1)
235    #     SRRF_array = pps[0:height_m, 0:width_m]
236
237    #     return SRRF_array
238
239
240    # def calculate_pairwise_product_sum_(self, rad_array):
241    #     n_time_points, height_m, width_m = rad_array.shape
242
243    #     SRRF_array = np.zeros((height_m, width_m), dtype=np.float32)
244    #     max_array = np.max(rad_array, axis=0)
245
246    #     for i in range(height_m):
247    #         for j in range(width_m):
248    #             r0 = np.maximum(max_array[i][j], 0)  # maximum value at position (i,j)
249    #             r1 = np.maximum(max_array, 0)  # maximum values of all pixels in the image
250    #             pps = np.sum(r0 * r1) / np.maximum((n_time_points * (n_time_points + 1)) / 2, 1)
251    #             SRRF_array[i][j] = pps
252
253    #     return SRRF_array
254
255
256
257
258
259
260
261
262

def calculate_SRRF_temporal_correlations( im: <built-in function array>, order: int = 1, do_integrate_lag_times: bool = 0): View Source

 5def calculate_SRRF_temporal_correlations(im: np.array, order: int = 1, do_integrate_lag_times: bool = 0): #these are for SRRF
 6
 7    im = np.array(im, dtype='float32')
 8    assert im.ndim == 3 and order <= 4
 9    
10    if order == 0: #TODO: check order number in NanoJ
11        out_array = np.amax(im, axis=0)
12    
13    elif order == 1:
14        out_array = np.mean(im, axis=0)
15    
16    elif order == -1:
17        out_array = calculate_pairwise_product_sum(im)
18
19    else: # order = 2 or order = 3 or order = 4
20        out_array = calculate_acrf_(im, order, do_integrate_lag_times)
21    
22    return out_array

def calculate_eSRRF_temporal_correlations(im: <built-in function array>, correlation: str): View Source

25def calculate_eSRRF_temporal_correlations(im: np.array, correlation: str):
26
27    im = np.array(im, dtype='float32')
28
29    if correlation == "AVG":
30        out_array = np.mean(im, axis=0)
31
32    elif correlation == "VAR":
33        out_array = np.var(im, axis=0)
34
35    elif correlation == "TAC2":
36        out_array = calculate_tac2(im)
37
38    else:
39        raise ValueError(f"Type of correlation must be AVG, VAR or TAC2")
40    
41    return out_array

def calculate_pairwise_product_sum(rad_array): View Source

45def calculate_pairwise_product_sum(rad_array):
46    n_time_points, height_m, width_m = rad_array.shape
47
48    out_array = np.zeros((height_m, width_m), dtype=np.float32)
49    #max_array = np.max(rad_array, axis=0)
50    counter = 0
51    pps = 0
52
53    for t0 in range(n_time_points):
54        r0 = np.maximum(rad_array[t0],0)
55        if np.any(r0) > 0:
56            for t1 in range(t0, n_time_points):
57                r1 = np.maximum(rad_array[t1],0)
58                pps += r0 * r1
59                counter += 1
60        else:
61            counter += n_time_points - t0
62    pps = pps/ max(counter,1)
63    out_array = pps
64
65    return out_array

def calculate_acrf_(rad_array, order, do_integrate_lag_times): View Source

 68def calculate_acrf_(rad_array, order, do_integrate_lag_times):
 69    im = rad_array.copy()
 70    n_time_points, height_m, width_m = im.shape
 71    mean = np.mean(im, axis=0)
 72
 73    out_array = np.zeros((height_m, width_m), dtype=np.float32)
 74
 75    abcd = np.zeros((height_m, width_m), dtype=np.float32)
 76    abc = np.zeros((height_m, width_m), dtype=np.float32)
 77    ab = np.zeros((height_m, width_m), dtype=np.float32)
 78    cd = np.zeros((height_m, width_m), dtype=np.float32)
 79    ac = np.zeros((height_m, width_m), dtype=np.float32)
 80    bd = np.zeros((height_m, width_m), dtype=np.float32)
 81    ad = np.zeros((height_m, width_m), dtype=np.float32)
 82    bc = np.zeros((height_m, width_m), dtype=np.float32)
 83
 84    if do_integrate_lag_times != 1:
 85        t = 0
 86        while (t < n_time_points - order):
 87            ab = ab + (im[t] - mean) * (im[t+1] - mean)
 88            if order == 3:
 89                abc = abc + (im[t] - mean) * (im[t+1] - mean) + (im[t+2] - mean)
 90            if order == 4:
 91                a = im[t] - mean
 92                b = im[t+1] - mean
 93                c = im[t+2] - mean
 94                d = im[t+3] - mean
 95                abcd = abcd + np.multiply(np.multiply(np.multiply(a,b),c),d)
 96                cd = cd + np.multiply(c,d)
 97                ac = ac + np.multiply(a,c)
 98                bd = bd + np.multiply(b,d)
 99                ad = ad + np.multiply(a,d)
100                bc = bc + np.multiply(b,c)
101            t = t + 1
102        if order == 3:
103            out_array = np.absolute(abc) / n_time_points
104        elif order == 4:
105            out_array = np.absolute(abcd - ab * cd - ac * bd - ad * bc) / n_time_points
106        else:
107            out_array = np.absolute(ab) / n_time_points
108    
109    else:
110        n_binned_time_points = n_time_points
111        tbin = 0
112        while n_binned_time_points > order:
113            t = 0
114            ab = np.zeros((height_m, width_m), dtype=np.float32)
115            while (t < n_binned_time_points - order):
116                tbin = t * order
117                ab = ab + (im[t] - mean) * (im[t+1] - mean)
118                if order == 3:
119                    abc = abc + (im[t] - mean) * (im[t+1] - mean) + (im[t+2] - mean)
120                if order == 4:
121                    a = im[t] - mean
122                    b = im[t+1] - mean
123                    c = im[t+2] - mean
124                    d = im[t+3] - mean
125                    abcd = abcd + np.multiply(np.multiply(np.multiply(a,b),c),d)
126                    cd = cd + np.multiply(c,d)
127                    ac = ac + np.multiply(a,c)
128                    bd = bd + np.multiply(b,d)
129                    ad = ad + np.multiply(a,d)
130                    bc = bc + np.multiply(b,c)
131                    
132                im[t] = np.zeros((height_m, width_m), dtype=np.float32)
133
134                if tbin < n_binned_time_points:
135                    for _t in range(order-1):
136                        im[t] = im[t] + np.divide(im[tbin + _t], order)
137                t = t + 1
138            if order == 3:
139                out_array = np.absolute(abc) / n_binned_time_points
140            elif order == 4:
141                out_array = np.absolute(abcd - ab * cd - ac * bd - ad * bc) / n_binned_time_points
142            else:
143                out_array = np.absolute(ab) / n_binned_time_points
144
145            n_binned_time_points = n_binned_time_points / order
146
147    return out_array

def calculate_tac2(rad_array): View Source

150def calculate_tac2(rad_array):
151    mean = np.mean(rad_array, axis=0)  
152    centered = rad_array - mean  # center data around the mean
153    nlag = 1  # number of lags to compute TAC2 for
154    out_array = np.mean(centered[:-nlag] * centered[nlag:], axis=0)
155
156    return out_array
157
158
159
160
161# class TemporalCorrelation:
162#     accepted_correlation_types = ('mean', 'var', 'tac2', 'mip', 'pair_wise_product_sum')
163
164#     def __init__(self, correlation_type: str):
165#         """
166#         Perform a temporal correlation analysis of an image (with shape time, height, width) #TODO: discuss what shape should be the input
167#         :param correlation_type: desired type of interpolation ("mean", "var" or "tac2")
168#         """
169#         self.correlation_type = correlation_type
170
171
172#     def calculate_tc(self, im: np.array):
173
174#         out_array = np.empty((im.shape[0]))
175
176#         # assert isinstance(im, np.ndarray)
177
178#         if im.dtype != np.float32:
179#             im = np.array(im,dtype='float32')
180
181#         assert im.ndim == 3 #TODO: discuss if we should consider (t,c,z,r,c) here
182
183#         if self.correlation_type == "mean":
184#             out_array = np.mean(im, axis=0)
185        
186#         elif self.correlation_type == "var":
187#             out_array = np.var(im, axis=0)
188        
189#         elif self.correlation_type == "tac2": # second order autocorrelation function
190#             mean = np.mean(im, axis=0)  
191#             centered = im - mean  # center data around the mean
192#             nlag = 1  # number of lags to compute TAC2 for
193#             out_array = np.mean(centered[:-nlag] * centered[nlag:], axis=0)
194
195#         elif self.correlation_type == "mip": # maximum intensity projection
196#              out_array = np.amax(im, axis=0) #calculate_mip(im, doIntensityWeighting)
197
198#         elif self.correlation_type == "pair_wise_product_sum":
199#             out_array = calculate_pairwise_product_sum_(im)
200        
201#         else:
202#              raise ValueError(f"Type of correlation must be one of {self.accepted_correlation_types}")
203
204#         return out_array
205    
206    
207    # def calculate_mip(im: np.array, doIntensityWeighting: bool = True):
208    #     """
209    #     Calculate Maximum Intensity Projection of an input image
210    #     """
211    #     #nb_time_points, height, width = im.shape
212    #     out = np.empty((im.shape[0]))
213
214    #     if doIntensityWeighting == True:
215    #         out = 
216    #     else:
217    #         out = np.amax(im, axis=0)
218    #     return out
219
220    # def calculate_pairwise_product_sum(im: np.array):
221    #     """
222    #     Calculate Pair-Wise Product Sum of a 3D input image
223    #     """
224
225
226
227    # def calculate_pairwise_product_sum(self, rad_array: np.array):
228    #     n_time_points, height_m, width_m = rad_array.shape
229
230    #     SRRF_array = np.zeros((height_m, width_m), dtype=np.float32)
231
232    #     r0 = np.maximum(rad_array[:, :, None, :], 0) 
233    #     r1 = np.maximum(rad_array[:, :, :, None], 0) 
234
235    #     pps = np.sum(r0 * r1, axis=(2, 3)) / np.maximum((n_time_points * (n_time_points + 1)) / 2, 1)
236    #     SRRF_array = pps[0:height_m, 0:width_m]
237
238    #     return SRRF_array
239
240
241    # def calculate_pairwise_product_sum_(self, rad_array):
242    #     n_time_points, height_m, width_m = rad_array.shape
243
244    #     SRRF_array = np.zeros((height_m, width_m), dtype=np.float32)
245    #     max_array = np.max(rad_array, axis=0)
246
247    #     for i in range(height_m):
248    #         for j in range(width_m):
249    #             r0 = np.maximum(max_array[i][j], 0)  # maximum value at position (i,j)
250    #             r1 = np.maximum(max_array, 0)  # maximum values of all pixels in the image
251    #             pps = np.sum(r0 * r1) / np.maximum((n_time_points * (n_time_points + 1)) / 2, 1)
252    #             SRRF_array[i][j] = pps
253
254    #     return SRRF_array